The wave packet propagation using wavelets 
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Abstract 

It is demonstrated that the wavelets can be used to considerably speed up simulations 
of the wave packet propagation in multiscale systems. Extremely high efficiency is 
obtained in the representation of both bound and continuum states. The new method 
is compared with the fast Fourier algorithm. Depending on ratios of typical scales of a 
quantum system in question, the wavelet method appears to be faster by a few orders 
of magnitude. 

1. Owing to the fast development of the computational tools the direct solution of the 
time-dependent Schroedinger equation has become one of the basic approaches to study the 
evolution of quantum systems. Thus, the wave packet propagation (WPP) method is success- 
fully applied to time dependent and time-independent problems in gas-surface interactions, 
molecular and atomic physics, and quantum chemistry ^ ^, One of the main issues of 
the numerical approaches to the time-dependent Schroedinger equation is the representation 
of the wave function \^{t)) of the system. Development of the pseudospectral global grid 
representation approaches yielded a very efficient way to tackle this problem. The discrete 
variable representation (DVR) and the Fourier grid Hamiltonian method (FGH) |^ |^ 
have been widely used in time-dependent molecular dynamics [|l], H, ||, S-matrix p], |10| or 
eigenvalue 0] calculations. The FGH method based on the Fast Fourier Transform (FFT) 
algorithm is especially popular. This is because for the mesh of points the evaluation of 
the kinetic energy operator requires only NlogN operations and can be easily implemented 
numerically. In the standard FGH method the wave function is represented on a grid of 
equidistant point in coordinate and momentum space. It is well suited when the the local 
de Broglie wavelength does not change much over the physical region spanned by the wave 
function of the system |]lT| . There are, however, lot of examples where the system under con- 



sideration spans the physical regions with very different properties. Consider, for example, 
the scattering of a slow particle on a narrow and deep potential well. Despite the short wave 
lengths occur only in the well, in the FGH method they would determine the lattice mesh 
over entire physical space leading to high computational cost. In fact, pseudospectral global 
grid representation approaches are difficult to use in multiscale problems. 

This is why much of work has been devoted recently to the development of the mapping 
procedures in order to enhance sampling efficiency in the regions of the rapid variation of 
the wave function |TT], |T^, |T^, |T^, 0. Though mapping procedure, based on the variable 



change x = /(^) is very efficient in ID, it is far from being universal. One obvious drawback 
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is that it is difficult to implement in higher dimensions. In the case of several variables, there 
is no simple procedure to define the mapping functions / so that the lattice would be fine only 
in some designated regions of space. Next, the topology of the new coordinate surfaces can 
be different from that of Cartesian planes. Therefore the Jacobian may vanish at some points 
(e.g., spherical coordinates). This leads to singularity in the kinetic energy and imposes small 
time step in simulations [|I^]. 

In this letter we propose a novel approach to the wave packet propagation in multiscale 
systems. It is based on the use of wavelets as basis functions in the projected Hilbert space. 
In contrast to the plane waves, the wavelets are localized in both real and momentum spaces. 
This characteristic property of wavelets allows us to accurately describe short wave length 
components of the wave function in designated spatial regions, while keeping only few basis 
elements in the bulk. This leads to a drastic reduction of the projected Hilbert space dimension 
without any loss of accuracy. In some cases it may even become possible to diagonalize the 
Hamiltonian matrix so that the time evolution can be elementary followed in the basis of 
eigenstates. We illustrate our approach by simulations of a simple multiscale one-dimensional 
scattering problem. A generalization to higher dimensions is straightforward and based on 
the standard mathematical construction of multidimensional wavelets [T^. Wavelet bases 



have been successfully used for the systematic treatment of the multiple length scales in 
the electronic structure of matter |18|, |19|, |2ll] (for a comprehensive review see and 



references therein). Here we apply, for the first time, wavelets to study the time evolution in 
multiscale quantum systems and demonstrate the advantages of the wavelet method over the 
conventional FGH method. 

2. To illustrate the efficiency of the wavelet method in the wave packet propagation, we 
consider the ID scattering of an electron on a narrow and deep potential well. Despite we 
use an electron as a projectile, procedure described below readily applies to molecular or 
atomic wave packets. The potential has the form V = — 64exp[— (lOOx)^] (atomic units are 
used throughout the paper) and supports single bound state with the energy —0.63200. We 
are interested in the transmission and refiection coefficients for the energies of the scattered 
electron within 0.2 — 1.1 range. This corresponds to a typical wave length A = 27r which is 
much larger than the potential well width. To calculate the scattering properties of the system 
we use the Gaussan wave packet impinging on a potential well. The calculation requires a 
typical size of the box x e [—60, 60], where 18 a.u. from each side are allocated for an optical 
potential Vopt that absorbs reflected and transmitted waves 0]. 

We flrst introduce a uniform lattice and simulate the wave packet propagation by means 
of the FGH method. Split-operator technique |]^, |TB[ is used to calculate the action of the 
evolution operator (W = V + Vopt)'- 



|^(t + At)) = ex.p{-iAt H + Vopt)\'^it)) 



exp{-iAtW/2) exp{-iAtf) exp{-iAtW /2) |^(t)) + 0{Af 



(1) 



Numerical convergence is obtained with = 2^^ points at the grid and a time step At = 
0.0002. Calculated transmission and reflection coefficients are shown in Fig. 1. The time 



2 



evolution of the wave packet is presented in Fig.2. The colors represent the magnitude of the 
wave packet in the logarithmic scale: —10 > In |\l/()f:, x)| > 0.3. As the color changes from red 
to violet, In decreases from 0.3 to —10. Without the potential present, we would see 

a colored ray (a trajectory of a free particle is a straight line) spreading for larger values of t. 

To compare the above results with our new wavelet method, we take the Daubechies 
wavelets Pio which have ten first vanishing moments and the filter of length 20 The 
basis is generated by the scaling function (j){x) whose support lies in the interval x G [0, 19]. 
The lowest resolution level is set so that the initial wave packet is reproduced with high 
accuracy in the orthonormal basis of the functions = ^/2(f)(2x — j) where j runs over 

integers inside the interval (—120,120]. The wavelets ifjkjix) = ^/^■\\}{2^x — j) are used to 
describe short wave length components of the wave function in the vicinity of the potential 
well. For a moment, they should be regarded as a special orthonormal basis with compact 
support and the following properties. If the basis of scaling functions 0i spans well 
functions whose Fourier transforms are concentrated in a wave length band around Ai, the 
corresponding wavelets from the kth resolution level form a good basis for a band centered 
at 2-'=Ai. In our case we have used k = 1,2,3,4. The number of wavelets needed on each 
resolution level is determined by the potential width (by its shape, in general). The technical 
details are explained in Section 3. It appears necessary to take 20 wavelets on each resolution 
level. Thus, all together we use only 240 + 4-20 = 320 coefficients (or, equally basis functions) 
to describe the wave packet propagation instead of 2^^ in the fast Fourier method. Due to such 
a tremendous reduction of the size of the problem, the Hamiltonian matrix can be directly 
diagonalized and the time dependent wave function can be easily obtained as: 

\^{t + At)) = e-'^°-'^' Ee e"'""^* 1^) {E\^{t)) , (2) 

where \E) and E stands for the eigenvector and eigenvalue of the Hamiltonian, respectively. 
Convergent results are obtained with time step dt = 0.025, i.e., much larger than in the fast 
Fourier method above. The results are given on Figs. 1 and 3. There is a perfect coincidence 
of the fast Fourier and our wavelet results, while the wavelet method is nearly 100 times 
faster. 

For the sake of comparison, we also used Eq. (1) to simulate the time evolution. The same 
accuracy is achieved for a time step comparable to the time step in FGH method. However 
there is an alternative, much better splitting scheme thanks to the localization properties of 
wavelets. Consider the decomposition H = H1 + H2 where Hi contains matrix elements of the 
basis elements localized in the vicinity of the potential well, while H2 corresponds to the free 
space and, therefore, contains matrix elements only of the scaling functions. The propagation 
can then be done accordingly: 

\^{t + At)) = e-*^°p*^*e-'^^^*/2 l^^^e-'^i^* \Ei) {Ei\] e''^'^'^^ \^{t)) , (3) 

where \Ei) and Ei stands for the eigenvector and eigenvalue of the Hi. In the spit method 
([^) the error depends on the spectral range of operators involved. For larger spectral ranges, 
the errors get larger. The conventional way to cope with this problem is to reduce the time 
step. The approach (^, which we call the "wavelet tower diagonalization" , offers a better 
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alternative. The convergence here is drastically improved because (a) [Hi,H2] is small and 
(b ) the operator Hi with a large spectral range can be diagonalized. Both the properties are 
hardly achievable without the wavelet basis. In our example the convergence is reached for 
a time step dt = 0.01 (vs dt = 0.0002 in ([I|)). Corresponding results are presented in Fig. 
1. This approach applies even better to evolving heavy-particle wave packets because the 
Chebyshev |T^, ^ or Lanczos [0, ^| schemes can be used. These schemes are, first, more 



efficient than the split method and, second, allow one to take full advantage of the sparse 
structure of the Hamiltonian in the wavelet basis (see below). 

Note also that by making the potential depth greater so that the minimal wave length 
becomes twice shorter, we would have to increase the lattice size by factor two in the Fourier 
method, thus increasing the number of coefficients needed to describe the wave packet by 2^^ 
(!), while in our wavelet approach one more resolution level is to be added, implying only 
20 extra wavelet coefficients in the wave packet decomposition. Thus, the wavelet approach 
offers a systematic and easy way to improve the accuracy of simulations. 

3. The description of technical details of our approach is limited to essential practical 
steps necessary to reproduce our results and to apply the algorithm to new systems. A 
general theoretical analysis of wavelets bases in multidimensional spaces can be found in |r7|. 



So we only discuss the Hilbert space L2 of square integrable functions of a real variable x, 
J dx\'i'{x)\'^ < 00. 

(i). The scaling function is defined by the equation (j){x) = 2Y^i hi(j){2x — I) (called 
the scaling relation), where the coefficients hi are called a filter. For a finite filter, the scaling 
function has compact support. Define (f>n,jix) = 2"/^0(2"'x — j) for all integers j and n. The 
filter hi satisfies an equation obtained by combining the scaling relation with the required 
orthonormality condition of the scaling functions 0„ j for fixed n. Given a filter, numerical 
values of (f){x) can be generated by an iteration procedure Jl^ ^ . 



(ii). The subspace of L2 spanned by 0„j is denoted Vn- An important property Vn is 
that Vn-i C Vn and the projection of any function \l/ from L2 onto Vn converges to \l/ in the 
L2 sense as n ^ 00. Consider an orthogonal decomposition: Vn+i = Vn® Wn- There exists 
a function G Wq called the mother wavelet such that -ipnjix) = 2"/^?/'(2"'x — j) form an 
orthonormal basis in Wn- The numerical values oiip can be generated from the scaling relation 
ip^x) = 2 X)z(— l)'^i-/0(2x — I). A finite dimensional subspace of Vn is used to approximate 
L2 in simulations. By construction, the functions (pij, I < n, and tpkj for k = 1,1 + 1, ...,n — 1 
form an orthonormal basis in K: {i!k,j\ipk',i) = hw^ij, {4>iJ(j)i,i) = 6ij, and {ipk,j\4>i,i) = 0. 

(in). From the definition of (pij and ipk,j it is clear that the index j indicates the position 
of support of the scaling function or wavelet. The index k of ip^j can be understood as follows. 
Let the Fourier transform of the mother wavelet ip be peaked at a momentum p^. Then the 
Fourier transform of ipk,j is peaked at the momentum = 2^pq. Since the wavelets with 
different k are orthogonal, the coefficients dkj = {iphjl"^) determine relative amplitudes of 
successively shorter wave length components of \l/ in the vicinity of x = 2~'^j as k increases. 
For this reason, the index j is called a position index, and k is called a resolution level. 

(iv). From the physical properties of a system one can estimate a wave length band, 
A G [Xmin, ^max], required for simulations. The lowest resolution level / is identified by the 
condition that (j)ij span functions with Fourier components in the vicinity of Xmax- The 
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necessary maximal resolution level is determined hj n ^ ^og2{Xmax / ^min) ■ In our example 
/ = 1 and n = 5. If scaling functions are needed (to cover a given physical volume), then, 
in general, on every wavelet resolution level there will be Jk = 2'^~^N basis functions. The 
wave packet is decomposed in the corresponding basis of Vn 

N n-1 Jk 

*W = E siAt)(t>i, + E E dk,{t)i'k, . (4) 
j=i k=i j=i 

The decomposition coefficients as functions of time are to be found from the Schroedinger 
equation. Yet, the total number of coefficients in is about the same as that in the uniform 
finite lattice approach. 

The advantage of using wavelet bases becomes significant whenever the volume of regions 
where short wave lengths can appear during the time evolution is much smaller than the total 
physical volume of the system. This allows one to substantially reduce in Eq.(^) by taking 
higher resolution level wavelets only where they are needed. In our case we need only 20 
wavelets at each resolution level. We shall refer to all the wavelets needed in the vicinity of 
one local minima of the potential as a wavelet tower. Higher tower fioors correspond to higher 
resolution levels k. Each wavelet is thought of as a building block of width equal . 

(v). An important part of the algorithm is the projection of the Hamiltonian onto the 
wavelet towers. Matrix elements of the potential as well as the initial coefficients sij = 
{(pijl"^), dkj = (iphjl"^) are computed via Riemann sums. The matrix elements of the kinetic 
energy are computed using the the fast Fourier transform to obtain the second derivative 
of the function. The Hamiltonian matrix is sparse because of compact support of the basis 
functions. For instance, in our example {(j)n,j\H\(f)n,i) = 0, if \i — j\ > 20. Finally, the time 
evolution can be computed by standard techniques such as split, Lanczos, or Chebychev 



methods ||T6| , |23| , p4| . If after the projection on wavelet towers, the Hamiltonian matrix is 
small enough for a direct diagonalization, Eq. (Q) can be used for the propagation. 

In conclusion, we have demonstrated that wavelet bases can be extremely efficient for 
solving the time- dependent Schroedinger equation for multiscale systems. The bound and 
continuum states are accurately described with a much less number of basis functions as 
would be required in the standard Fourier-grid methods. In our example we were able to 
reduce the basis size by a factor of 50, which allowed us to scale down the computation time 
by a factor of 100. The method can easily be implemented in higher dimensions where a 
wavelet basis is built by taking the direct product of one-dimensional wavelets |T^. The 
wavelet approach becomes more efficient as ratios of typical scales get larger. Finally we 
would like to stress that the wavelet method is ideally suitable for simulating wave packet 
propagation in multiscale problems with complicated form of the potential. This is because 
the wavelet towers can be custom designed for any topology of the potential minima (e.g., 
potential valleys) regardless of the dimensionality of the system. 
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Figure captions 

Fig. 1. Calculated transmission (black) and reflection (red) coefficients. Circles: The stan- 
dard approach based on the split propagation and Fourier grid with uniform mesh. Solid 
curves: The wavelet approach with Hamiltonian matrix diagonaliztion. Triangles: The split 
propagation with the "wavelet tower diagonalization" as in Eq. @. 
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